Flowpipe construction consists of under or over-approximating the sets of states reachable by dynamical systems. Recently a method has been developed for the class of polynomial ODEs with uncertain initial states (see [1], abbreviated XFZ18under). This method consists of reducing the Hamilton-Jacobi-Bellman equation to a hierarchy of semidefinite programs.
In this notebook we consider the problem of approximating the flowpipe of a system of polynomial ODEs using XFZ18under. This is a Julia implementation that relies on the JuMP ecosystem (JuMP, PolyJuMP, SumOfSquares, MathProgInterface, MathOptInterfaceMosek) and the JuliaAlgebra ecosystem (MultivariatePolynomials, DynamicPolynomials). The implementation is evaluated on a set of standard benchmarks from formal verification and control engineering domains.
References:
Quick links to documentation:
Consider the following van-der-Pol system:
$$ \dot{x}_1 = x_2 \\ \dot{x}_2 = -0.2x_1 + x_2 - 0.2x_1^2 x_2 $$
using MultivariatePolynomials,
JuMP,
PolyJuMP,
SumOfSquares,
DynamicPolynomials,
MathOptInterfaceMosek,
MathematicalSystems
const ∂ = differentiate
# symbolic variables
@polyvar x₁ x₂ t
# time duration (scaled, see dynamics below)
T = 1.0
# dynamics
f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂]
# set of initial states X₀ = {x: V₀(x) <= 0}
V₀ = x₁^2 + x₂^2 - 0.25
# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]
g = 9 - x₁^2 - x₂^2
# degree of the relaxation
k = 4
# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i
X = monomials([x₁, x₂], 0:k)
XT = monomials([x₁, x₂, t], 0:k)
# create a SOS JuMP model to solve with Mosek
model = SOSModel(with_optimizer(MosekOptimizer))
# add unknown Φ to the model
@variable(model, Φ, Poly(XT))
# jacobian
∂t = α -> ∂(α, t)
∂xf = α -> ∂(α, x₁) * f[1] + ∂(α, x₂) * f[2]
LΦ = ∂t(Φ) + ∂xf(Φ)
# Φ(x, t) at time 0
Φ₀ = subs(Φ, t => 0.)
# scalar variable
@variable(model, ϵ)
@variable(model, s₁, SOSPoly(XT))
@variable(model, s₂, SOSPoly(XT))
@variable(model, s₄, SOSPoly(XT))
@variable(model, s₅, SOSPoly(XT))
@variable(model, s₇, SOSPoly(X))
@variable(model, s₉, SOSPoly(X))
@constraint(model, ϵ >= 0.)
@constraint(model, LΦ - s₁*t*(T-t) - s₂*g ∈ SOSCone())
@constraint(model, ϵ - LΦ - s₄*t*(T-t) - s₅*g ∈ SOSCone())
@constraint(model, Φ₀ - V₀ - s₇*g ∈ SOSCone())
@constraint(model, ϵ + V₀ - Φ₀ - s₉*g ∈ SOSCone())
@objective(model, Min, ϵ)
optimize!(model)
println("Relaxation order : k = $k")
println("JuMP.termination_status(model) = ", JuMP.termination_status(model))
println("JuMP.primal_status(model) = ", JuMP.primal_status(model))
println("JuMP.dual_status(model) = ", JuMP.dual_status(model))
println("JuMP.objective_bound(model) = ", JuMP.objective_bound(model))
println("JuMP.objective_value(model) = ", JuMP.objective_value(model))
# Recovering the solution:
ϵopt = JuMP.objective_value(model)
# Punder <= 0
Punder = subs(JuMP.value(model[:Φ]), t => T)
# Pover <= 0
Pover = subs(JuMP.value(model[:Φ]), t => T) - ϵopt * (T+1)
using ImplicitEquations, Plots
gr() # or pyplot()
G = plot()
_Punder(x, y) = sum([Punder.a[i]*x^exponents(mi)[1]* y^exponents(mi)[2] for (i, mi) in enumerate(Punder.x)])
plot!(G, _Punder ⩵ 0., xlims=(-10, 10), ylims=(-10, 10), color="red")
_Pover(x, y) = sum([Pover.a[i]*x^exponents(mi)[1] * y^exponents(mi)[2] for (i, mi) in enumerate(Pover.x)])
plot!(G, _Pover ⩵ 0., xlims=(-10, 10), ylims=(-10, 10), color="blue")
G
G = plot()
_Punder(x, y) = sum([Punder.a[i]*x^exponents(mi)[1]* y^exponents(mi)[2] for (i, mi) in enumerate(Punder.x)])
Gu = plot(_Punder ≪ 0., xlims=(-8, 8), ylims=(-8, 8))
_Pover(x, y) = sum([Pover.a[i]*x^exponents(mi)[1] * y^exponents(mi)[2] for (i, mi) in enumerate(Pover.x)])
Go = plot(_Pover ≪ 0, xlims=(-8, 8), ylims=(-8, 8))
plot(Gu, Go)
| Package | k | Constraints | Scalar variables | Matrix variables | Time(s) |
|---|---|---|---|---|---|
| JuMP | 2 | 1017 | 803 | 10 | < 1 |
| YALMIP | 2 | 152 | 63 | 10 | < 1 |
| JuMP | 3 | 2871 | 2471 | 10 | ~ 1 |
| YALMIP | 3 | 254 | 121 | 10 | ~ 1 |
| JuMP | 4 | 7119 | 6450 | 10 | 6.22 |
| YALMIP | 4 | 394 | 206 | 10 | 1.18 |